
clear
do "...\First.do"

********************************************************************************
use "$work\analysis_sample.dta"
keep if inrange(age,40,70)

********************************************************************************

** Average deaths per relative time by group 
bys timing: egen mean_ll=mean(death) if low_ses==1 & help3==1
bys timing: egen mean_lh=mean(death) if low_ses==1 & help3==0
bys timing: egen mean_hl=mean(death) if low_ses==0 & help3==1
bys timing: egen mean_hh=mean(death) if low_ses==0 & help3==0

bys low_ses help3: g n=_N/10

g upperll=mean_ll2+1.96*sd_ll2/sqrt(n)
g lowerll=mean_ll2-1.96*sd_ll2/sqrt(n)

g upperlh=mean_lh2+1.96*sd_lh2/sqrt(n)
g lowerlh=mean_lh2-1.96*sd_lh2/sqrt(n)

g upperhl=mean_hl2+1.96*sd_hl2/sqrt(n)
g lowerhl=mean_hl2-1.96*sd_hl2/sqrt(n)

g upperhh=mean_hh2+1.96*sd_hh2/sqrt(n)
g lowerhh=mean_hh2-1.96*sd_hh2/sqrt(n)


preserve 
keep timing mean_hh mean_hl mean_lh mean_ll mean_hh2 mean_hl2 mean_lh2 mean_ll2 upper* lower*
duplicates drop

foreach var of varlist mean_hh mean_hl mean_lh mean_ll mean_hh2 mean_hl2 mean_lh2 mean_ll2 upperll lowerll lowerlh upperlh upperhl lowerhl upperhh lowerhh {
    replace `var'=`var'*100
}

set scheme s2mono
sort timing

tw (line mean_hh timing, lcolor(sienna)) (line mean_hl timing, lcolor(navy)) ///
(line mean_ll timing, lcolor(sienna)) (line mean_lh timing, lcolor(navy)  graphregion(color(white)) bgcolor(white) legend(order(1 "Patient high-SES - Physician high-SES" 2 "Patient high-SES - Physician low-SES" 3 "Patient low-SES - Physician low-SES" 4 "Patient low-SES - Physician high-SES") region(lwidth(none)) col(1) pos(11) ring(0) size(small)) ytitle("Mortality %") xtitle("Years since clinic closure") xscale(range(-4 5)) xlabel(-4(1)5) xline(0, lp(dash) lc(*0.5)))
graph export "$fig\Event_raw.png", replace 
graph export "$fig\Event_raw.pdf", replace 

restore

